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Abstract 

We present an equation-free multi-scale approach to the computational study of the collective 
dynamics of the Kuramoto model [Chemical Oscillations, Waves, and Turbulence, Springer- Verlag 
(1984)], a prototype model for coupled oscillator populations. Our study takes place in a reduced 
phase space of coarse-grained "observables" of the system: the first few moments of the oscillator 
phase angle distribution. We circumvent the derivation of explicit dynamical equations (approximately) 
governing the evolution of these coarse-grained macroscopic variables; instead we use the equation- 
free framework [Kevrekidis et al, Comm. Math. Sci. 1(4), 715 (2003)] to computationally solve 
these equations without obtaining them in closed form. In this approach, the numerical tasks for the 
conceptually existing but unavailable coarse-grained equations are implemented through short bursts 
of appropriately initialized simulations of the "fine-scale", detailed coupled oscillator model. Coarse 
projective integration and coarse fixed point computations are illustrated. 
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1 Introduction 



Coupled nonlinear oscillators can exhibit spontaneous emergence of order, a fundamental 
qualitative feature of many complex dynamical systems [Manrubia et al, 2004]. The col- 
lective, coarse-grained dynamics of coupled oscillator populations can provide insights in 
synchronization phenomena in many biological, chemical, and physical systems [Winfree, 
1967; Walker, 1969; Buck, 1988; Gray et al, 1989; Ertl, 1991; Wiesenfeld et al, 1996; 
Neda et al, 2000; Oliva & Strogatz, 2001; Kiss et al, 2002]. In certain prototype mod- 
els some exact results for the collective dynamics, mainly for an order parameter at the 
asymptotic states, can be obtained for either a few or an infinite number of oscillators 
(the so-called continuum limit) [Kuramoto, 1975; Kuramoto, 1984; Ariaratnam & Stro- 
gatz, 2001]; for instance, the continuum-limit Kuramoto model is well known to exhibit a 
phase transition at a critical coupling strength. However, general coupled oscillator models 
are seldom amenable to mathematical analysis, even for simple coupling topologies, and 
the understanding of the coarse-grained (emergent or macroscopic) dynamics of a given os- 
cillator population remains limited. Even fundamental questions, such as global stability, 
spectral properties of the stationary solution, and finite-number fluctuations, often remain 
unanswered [Strogatz, 2000]. 

Understanding the coarse-grained behavior of dynamical systems consisting of a large 
number of mutually interacting entities is one of the main goals of statistical mechanics. 
In some cases, a fine-scale (or microscopic) model with a large number of degrees of free- 
dom can be systematically reduced to a low-dimensional, coarse-grained model describing 
the dynamics by a small number of macroscopic variables (observables) . However, such a 
reduction is rarely available in contemporary science or engineering dealing with specific 
complex systems. One often has to deal with a situation where a reliable model is available 
only at a fine-scale level (atomistic, individual-based), yet he/she is practically interested in 
the emergent, system level, coarse-grained dynamics. Even if a useful coarse-grained model 
is believed to exist at a practical closure level, explicit formulas for this model are often 
unavailable; the systematic derivation of closed formulas for such models is difficult with- 
out introducing many simplifying assumptions. In some cases, it is also possible that the 
description of a system simply does not close at a particular level of coarse-graining. In the 
case that a coarse-grained description exists, yet it is unavailable in closed form, the recently 
developed equation-free framework provides a systematic, computer-assisted approach to 
exploring coarse-grained macroscopic behavior without explicit knowledge of the coarse- 
grained equations themselves [Theodoropoulos et al, 2000; Gear et al, 2002; Kevrekidis 
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et ai, 2003; Kevrekidis et ai, 2004]. This approach has been successfully combined with 
several classes of microscopic models, such as lattice-Boltzmann, kinetic Monte Carlo, and 
molecular dynamics simulators, to study the coarse-grained behavior of heterogeneous cat- 
alytic reactions, liquid crystal rheology, peptide fragment folding, etc. [Theodoropoulos et 
ai, 2000; Gear et ai, 2002; Makeev et ai, 2002; Siettos et ai, 2003, Hummer & Kevrekidis, 
2003]. In this paper we illustrate the equation- free approach to the computational explo- 
ration of the coarse-grained dynamics of the Kuramoto model at the strong coupling limit. 

The rest of the paper is organized as follows: The model under consideration and the 
regime in which we study it are described in Sec. 121 Some observations about the coarse- 
grained behavior and the equation-free analysis of phase-locked oscillator dynamics is pre- 
sented in Sec. |3] Conceptually, the generalization of this approach for different oscillator 
models is straightforward; possible shortcomings and limitations will be addressed in Sec. 
along with future research directions. 

2 Background; the Kuramoto model 

The Kuramoto model [Kuramoto, 1975] consists of N equally weighted, all-to-all, phase- 
coupled limit-cycle oscillators, where each oscillator has its own natural frequency uji drawn 
from a prescribed distribution function g(uj). We choose g{oj) to be a Gaussian distribution 
of standard deviation 0.1; however, our analysis is not limited to this particular choice. The 
microscopic individual level dynamics is easily visualized by imagining oscillators as points 
running around on the unit circle. Due to rotational symmetry, the average frequency 
Vt = YsiLi^i/N can be set to without loss of generality; this corresponds to observing 
dynamics in the co-rotating frame at frequency Q. The governing equation for the ith 
oscillator phase angle 9i is given by 

d0- K N 

^ = »i + x y LM0j-0i), i<*<n, (i) 

where K > is the coupling strength; in this paper we arbitrarily set N = 300 unless 
stated otherwise. 

It is known that as K is increased from above some critical value K c , more and 
more oscillators start to get synchronized (or phase-locked) until all the oscillators get fully 
synchronized at another critical value of K tp (Fig. [IJ. In our choice of Q = 0, the fully 
synchronized state corresponds to an exact steady state of the "detailed" , fine-scale problem 
in the co-rotating frame. Such synchronization dynamics can be conveniently summarized 
by considering the fraction of the synchronized (phase-locked) oscillators as in Fig. Q and 
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conventionally described by a complex order parameter [Kuramoto, 1975] 



1 N 

# (2) 



where the radius r measures the phase coherence, and ip is the average phase angle. Some 
exact results on the synchronization dynamics in the continuum-limit, in terms of this order 
parameter, are available [Kuramoto, 1984; Strogatz, 2000]. 

The transition from full to partial synchronization: We begin by restating cer- 
tain facts about the nature of the second transition mentioned above, a transition between 
the full and the partial synchronization regime at K = K tp , in the direction of decreasing 
K. A fully synchronized state (a "detailed" steady state) in the continuum limit corre- 
sponds to the solution to the following equations, the mean field type alternate form of 
Eq. (HJ) [Kuramoto, 1975]: 

d9- 

—i = n = oj i + rKsm(il;-e i ), l<i<N, (3) 

(JjL 

where Q is the common angular velocity of the fully synchronized oscillators (which is set 
to in our case). It can be rewritten as 

^— ^ = shn> -#,,), l<2<iV, (4) 

rK 

where the absolute value of the right hand side (RHS) is bounded by unity. As K approaches 
K tp from above, the LHS for the "extreme" oscillator (the oscillator in a particular family 
that has the maximum value of |0 — first exceeds unity, and a real- valued solution 
to Eq. (jlj) ceases to exist. We obtain the detailed steady state solution to Eq. (JTJ) as 
a function of K for 300 oscillators by using AUTO2000 (see Ref. [37]), which is based 
on fixed point computation and pseudo-arclength continuation. The bifurcation diagram 
plotted in terms of the extreme oscillator phase angle reveals that the transition to the 
partially synchronized state corresponds to a turning point bifurcation, of the so-called 
saddle node infinite period, or "sniper" type (Fig. EJ). Ermentrout and Kopell [Ermentrout 
and Kopell, 1984] have studied the transitions in the partial synchronization regime of a 
small number of oscillators, characterizing them as jumps in "frequency plateaus". 

Different random draws of u^s from g(u) for a finite number of oscillators result in 
slightly different values of K tp . We observe that K tp appears to follow a Gumbel type 
extreme distribution function [Kotz & Nadarajah, 2000], just as the maximum values of 
\Q — Ui\ do: 

p{K tp ) = a -l e -(Kt P -ti/* exp [_ e -(*«p-A0/<T] ; (5) 

where a and fi are parameters. 
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3 Coarse-grained dynamics of distribution moments 
3.1 Observation 

Our goal is to study the coarse-grained dynamics of an oscillator population, including 
the transient dynamics, in terms of coarse-grained observables. Since we are dealing with 
distributions, we naturally consider the first few moments as candidate coarse-grained vari- 
ables (i.e. a kinetic theory-like description). The moments of a distribution are often easily 
accessible in practice, and have a clear physical meaning. We consider the dynamics of the 
moments of the phase angle distribution function f(6), 



where n is a positive integer, and ( ) means an ensemble average. Since we choose to 
work with symmetric distributions, all the odd order moments are negligibly small. To 
obtain insights on the dynamics in this phase space, initial configurations of phase angles 
consistent with four sets of different values of the first two nonvanishing moments {M.2 and 
M.4) are generated. These initial phase angle configurations do not have any statistical 
correlations with the natural frequencies of the oscillators to which they are assigned; this 
is an important issue to which we will return below. Evolution of such configurations, 
obtained by direct numerical integration of Eq. reveals that the approach to a steady 
state consists of two stages. During the initial, relatively fast stage both moments decrease, 
and during a later slow stage the trajectory gradually approaches to an ultimate steady state 
along a common path, essentially a "slow manifold" of this dynamical system (Fig. |3J). The 
phase angle distributions after the initial stage, when the one-dimensional slow manifold 
finally has been approached, are very close to Gaussian distributions (marked by the dotted 
line in Fig. EJ). Once the trajectory approaches this slow manifold, the dynamics can be 
effectively described by (or observed on) any one even order moment. 

In the following (Sec. 13. 3J) . we will investigate the dynamics in the moment phase space 
following the equation- free approach [Theodoropoulos et ai, 2000; Kevrekidis et ai, 2003; 
Kevrekidis et ai, 2004]; this should be contrasted with the conventional approach of first 
attempting to derive closed governing equations for the leading moments, and subsequently 
analyzing these equations. 

3.2 Equation-free, coarse-grained multi-scale approach 

The basic idea of the equation-free approach is to utilize fine-scale simulation as an ex- 
periment, in the form of a time-stepper: short bursts of appropriately initialized fine-scale 





■5 



simulations are used to estimate quantities relevant for scientific computation with the 
coarse-grained model. Function evaluations or derivative evaluations for the coarse vari- 
ables involved in numerical tasks are thus substituted with estimation of these quantities, 
based on appropriately designed computational experiments with the detailed model. As 
a result, traditional continuum numerical techniques are directly "wrapped" around the 
fine-scale simulator. The essential steps can be summarized as follows: 

1. Identify a set of coarse-grained variables (observables) that sufficiently describe the 
system dynamics. We expect that a practically useful model exists and closes at the 
level of these coarse observables. In our analysis, we choose them to be the first 
two nonvanishing moments of the phase angle distributions, M2 and M.^. We will 
denote the microscopic description by u (u/s and 0's in our case) and the macroscopic 
description by U (.Mi's). 

2. Choose an appropriate lifting operator, /i, which maps the macroscopic description 
U to one or more consistent microscopic descriptions u; note that this mapping is not 
uniquely defined. This operator constructs one or more fine-scale realizations of the 
problem consistent with the prescribed values of the coarse observables U (u = /uU). 

3. Using a lifted, fine-scale realization as a new initial condition u(to), run the micro- 
scopic simulator to obtain the values u(t + AT) at a later time (AT > 0). This 
procedure may be repeated for many different realizations Uj(t ) consistent with the 
same macroscopic initial condition XJ(to), for purposes of ensemble averaging and vari- 
ance reduction. 

4. Use an appropriate restriction operator M which maps the microscopic description to 
the macroscopic description U(t) = M (u (t)). In our problem the restriction operator 
consists of simple evaluation of (ensemble averages of) the first two nonvanishing 
moments of the phase angle distribution, following their definitions in Eq. (jUJ). 

The above procedure results in the coarse timestepper <3>at: 

U(t„ + AT) = $ AT (U(t )) = M(u(to + AT)), (7) 

a map from coarse variables at time t — to to coarse variables at a later time t = to + AT, 
constructed through the microscopic simulation. Given such information, obtained for 
a number of appropriately chosen coarse initial conditions, it is possible to implement 
many scientific computing algorithms to solve the unavailable coarse-grained equations; 
see Sec. 13.31 for details of such computations using two different methods. An important 
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element of these algorithms is the fact that they function as protocols for the design of 
experiments with the fine-scale simulator: in the process of their execution they suggest 
new computational experiments with the coarse timestepper, and "bootstrap" their output 
toward the ultimate desired result. The success of this approach relies (as most of tradi- 
tional continuum numerical analysis) on smoothness of the dynamics with respect to the 
coarse variables of interest (in time, in space and in phase space) and on the identification 
of proper lifting and restriction operators. An extensive discussion of this method can be 
found in [Kevrekidis et at, 2003]. 

During the lifting step, we create phase angle distributions conditioned on a finite number 
of prescribed moments. This is accomplished as follows: Given g(ui), phase angles consistent 
with different values of the moments are generated by drawing random values from the 
inverse cumulative distribution of the following distribution function 

n(>2) 

1 + £ a p S p (9 2 ) , (8) 

p=2 

of which residual deviation from the Gaussian is expanded using Sonine polynomials (as- 
sociated Laguerre polynomials) [Chapman, 1970] as basis functions; f g {9 2 ) is a Gaussian 
distribution function with a prescribed value of M.2, and a p is the coefficient of the pth or- 
der Sonine polynomial S p . The value of M.2 P is determined only by a p , because the Sonine 
polynomials are otrhogonal with the Gaussian as the weighting function: 



fi(0) = W 2 ) 



oo 9 

e~ 9 S p {9 2 )S q {d 2 )d9 = 8 pq N p , (9) 



where 9 is generally a <i-dimensional vector, 5 pq is the Kronecker delta, and Af p is a 
normalization constant. The first three Sonine polynomials are given by Sq{9 2 ) = 1, 
S^O 2 ) = \d- 9 2 , and S 2 {9 2 ) = \d(d + 2) - \{d + 2)9 2 + \9\ where d is the dimen- 
sion, which in our case is unity. Our lifting operator takes us from the values of a few 
lower order moments to the values of an equal number of Sonine polynomial coefficients, 
and then to a consistent phase angle distribution through fi{9). We consider only even 
functions for fi{9), due to the symmetry of f(9) with our choice of g(u). 

One can clearly see in Fig. El that observation of the full system in the M.2 — -Ma phase 
space gives trajectories crossing one another; this suggests that a deterministic description 
closing with only these two variables is not possible. However, once the initial fast evolution 
stage has run its course, trajectories evolve on an apparently one-dimensional slow manifold, 
which suggests that the long-time coarse-grained dynamics could be described in terms of 
a single observable, say M.2 (or equivalently any other nonvanishing order moment). In 
other words, after sufficient time has elapsed, all the higher moments become slaved to 
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Ai 2 . Considering this observation, a more successful lifting operator may be constructed 
from the following two steps, when hA 2 is chosen to be a coarse variable of interest: (1) 
Drawing phase angles from an fi(0) satisfying a desired value Ai 2 = Ci, and (2) numerically 
integrating Eq. together with an algebraic constraint A4 2 — C\ = until the trajectory 
approaches the slow manifold. After this preparatory computation, the lifting step is 
complete; the constraint is then removed, and the system is allowed to evolve in order to 
evaluate its coarse timestepper. This lifting operation results in an augmented problem 
(ODEs with an algebraic constraint) that is described by a system of differential algebraic 
equations (DAEs) of index 1. For our illustrative example, it is straightforward to solve 
this set of DAEs by using the Lagrange multiplier method and the package DASSL (see 
Ref. [38 )• Direct integration of DAEs showing the relatively fast direct approach to the 
"mature" phase angle distributions on the slow manifold can be seen in Fig. Once the 
constrained trajectory approaches the slow manifold, the constraint is relaxed, and we can 
observe its subsequent approach to the steady state along the manifold. 

3.3 Equation-free computational results 

We now compute the steady state values of the moments in the fully synchronized regime 
(K > K tp ) using two different continuum numerical methods in the equation-free frame- 
work: coarse projective integration and coarse fixed point computation. 

We first evolve the system toward the stable steady state value of M. 2 using coarse 
projective integration [Gear & Kevrekidis, 2003a; 2003b]. The main assumption is that a 
macroscopic equation exists and closes in terms of Ai 2 ] our observations of the long-term 
dynamics are consistent with such an assumption. Starting from a relatively short time 
(At = 0.8) of direct integration, restriction of the results [shown as first group of five 
dots in Fig. |5] (a)] is used to estimate the coarse variable time derivative. Smoothness of 
the trajectory in time (Taylor series, as they appear in the simplest numerical integration 
schemes, such as forward Euler) is then used to "project" the value of the observable Ai 2 in 
forward time. We then lift with this estimated value [shown as the fist dot in the next group 
in Fig. (a)], and the whole procedure is repeated until the steady state is approached. 
During the evolution step, we do solve the full system of ODEs; but we do not solve it 
for all time: short bursts of detailed simulation are used to "jump" (and thus save) time, 
compared to a full integration [a solid line in Fig. |5] (a)]. 

As we discussed above, it is only after a relatively long time that the dynamics appear 
to "close" with only one observable Ai 2 ; this is because a certain evolution time is required 
before all higher moments of the distribution become slaved to Ai 2 . If we are interested in 
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iteration 


M 2 


\M 2 -Mf \/Mf 





l.OOOOxHT^ 


5.6768X10- 1 


1 


3.0362 xl0~ 2 


3.1261X10" 1 


2 


2.3179xl0~ 2 


2.0751 xl0~ 3 


3 


2.3096xl0- 2 


1.5131xl0" 3 



Table 1: Coarse fixed point computation at K — 0.7, using the Newton-Raphson method for M.2 only. 
The steady state value AA?, 3 , the variance of detailed steady state distribution, is 2.3131 x 10 -2 . The full 
system, Eq. was integrated for At = 2.0 at each iteration step, and nearby conditions separated by 
e = 1.0 x 10 -5 were used for linearization. 

dynamics occurring over shorter time scales, we need to work with a larger set of observ- 
ables. For example, if we are interested in dynamics over shorter times than those required 
for A^4 to get slaved to M.2-, we need to work with a coarse-grained model in terms of both 
M.2 and A^4 as coarse variables. The underlying premise is that the sequence of moments 
constitutes a singularly perturbed hierarchy; as longer times elapse, higher moments get 
progressively slaved to lower ones. Slow manifolds in moment space (graphs of functions, 
expressing higher moments as functions of lower ones) embody the closures that allow us 
to reduce the dimensionality of long-term dynamics. In equation-free computation, short 
simulation of the system itself is used to bring the system close to these slow manifolds (i.e., 
establish the closure numerically). At this level of coarse description, projective integration 
can again be done by lifting consistently with prescribed values of both M. 2 and A^4 (two 
constraints in the lifting DAE step) [Fig. El (b)]. 

Depending on the computational cost of projection, lifting, and restriction, these meth- 
ods have the potential for considerable computational savings in problems characterized by 
a large separation of time scales (gaps in their eigenvalue spectrum). In our case, we used 
a simple least squares fit algorithm for the estimation and subsequent projection in time, 
which required negligible computational effort. Restriction also required negligible compu- 
tation. However, each lifting step (integrating DAEs for times long enough to approach the 
slow manifold) could be relatively expensive. Working over shorter time scales (more coarse 
observables, higher dimensional slow manifolds) alleviates the duration of this preparatory 
lifting step. Finally, it is interesting to note that, in certain cases, forward integration 
can be used to evolve the coarse-grained system backward in time, possibly converging to 
certain types of saddle unstable states [Gear & Kevrekidis, 2004a]. 

In addition to coarse projective integration as a means of approaching the steady state 
in time, we also approximate it using a coarse contraction mapping, based on the Newton- 
Raphson method [Kevrekidis et ai, 2003]; we compute the solutions M^ ^ to the fol- 
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i.8000xi(r 2 


9.0000xl0~ 4 


5.8386x10-^ 
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2.3314xlCT 2 
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4.5518xl0~ 4 
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2.2902xl0~ 2 


1.6057xl0~ 3 


2.3052xl0" 4 



Table 2: Coarse fixed point computation at K = 0.7, using the Newton-Raphson method both for M 2 and 
Ma- The steady state value M% s , obtained from the detailed steady state distribution, is 1.6076 x 10~ 3 . 
The full system was integrated for At = 2.0 at each iteration step, and nearby conditions separated by 
e = 1.0 x 10~ 5 for M 2 (1.0 x 10 -6 for Mi) were used for linearization. 



lowing equations (of which explicit forms are unavailable): 

§at{M s 2 s[CG) )-M s 2 s{CG) = (10) 

or 

(CG) \ (CG) \ 



* AT ' Mf ca) ) ~ ( MT CG> I = °' (U) 



depending on the level of closure being considered, where the superscript stands for coarse- 
grained steady state. In this computation, residual and derivative evaluations required for 
each Newton-Raphson iteration are obtained by short bursts of direct integration of the 
microscopic model for AT = 2.0, starting at nearby initial conditions. Within only few 
Newton-Raphson iterations, an accurate solution to the unavailable equation can be found 
( Tables H] and |2I) ■ As a byproduct of this computation, upon convergence, the eigenvalues of 
the coarse linearization at steady state can be computed. In this case, when a steady state 
of the coarse problem corresponds to a true steady state of the fine-scale one, the eigen- 
values of the two are related (one expects to find, from the coarse procedure, the leading 
slow eigenvalues of the detailed problem linearization, discounting the neutral eigenvalue 
at that corresponds to rotational symmetry). As the number of coarse-grained variables 
increases, estimating each element of the coarse Jacobian becomes increasingly cumber- 
some; matrix-free methods of iterative linear algebra, like Newton-GMRES, are then used 
to solve the equations arising in the equation-free fixed point computation. These methods 
(such as Newton-GMRES; see the monograph [Kelley, 1995]) are based on the estimation of 
matrix-vector products, using short bursts of simulation with appropriately chosen pertur- 
bations of a given coarse initial condition (see also the Recursive Projection Method [Shroff 
& Keller, 1993]). 

4 Conclusions 

The long-term collective dynamic behavior of the Kuramoto model at the continuum (or 
thermodynamic) limit is often described in terms of an order parameter r [Strogatz, 2000] . 
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In this work we chose to observe the (long-time) collective dynamics for finite assemblies of 
coupled oscillators using a different set of coarse-grained variables: low order moments of 
the phase angle distribution. We found that, after sufficient evolution time has elapsed, the 
dynamics lie on a slow manifold parametrized by a few of these moments. The longer the 
dynamics have evolved, the lower the dimension of this manifold becomes, as higher order 
moments get progressively slaved to lower order ones. As a result, if we are interested 
in relatively long-term dynamics, we can work with less macroscopic observables in our 
lifting and restriction schemes. Two levels of closure were explored (parametrized by only 
one moment, and by two moments of the phase angle distribution respectively); shorter 
lifting computations were required for the two-moment description. We circumvented the 
derivation of closed governing equations for these coarse-grained variables, following the 
equation-free, coarse-grained multi-scale approach. Using this approach, we demonstrated 
coarse projective integration as well as coarse fixed point computation (followed by coarse 
stability analysis) in the fully synchronized regime. 

Our approach can, in principle, be extended in a straightforward manner to explore 
the dynamics of other coarse-grained variables of this model, or more complicated cou- 
pled oscillator models, as long as the long-term coarse-grained dynamics exhibit smooth 
low-dimensional behavior (i.e. are characterized by a slow, attracting manifold in the cor- 
responding phase space). A very important step is the identification of appropriate lifting 
and restriction operators, that allow us to (approximately) initialize the fine-scale system 
close to the slow manifold, and to compute the evolution of its macroscopic observables. 

As we saw in Fig. simply prescribing two moments of the distribution and using 
Eq. (jHJ) does not provide a good enough lifting; temporal trajectories of the coarse variables 
cross one another. We obtained a better lifting by introducing an additional preparatory 
step: evolution of the system constrained on the observables, until it approached the low- 
dimensional, slow manifold. This required the integration of a system of DAEs (in the 
spirit of algorithms like SHAKE [Ryckaert et ai, 1977] and umbrella sampling [Torrie & 
Valleau, 1974]). For this problem it was easy to augment the fine-scale model to introduce 
such a constraint; algorithms guiding to a point on the manifold without the imposition of 
an explicit constraint have also been devised, which was used for the detailed simulator as 
an unmodified "legacy code" [Gear & Kevrekidis, 2004b; Gear et ai, 2004c]. 

Suppose that the steady state phase angle distribution is known (e.g. from a long 
simulation), and we lift consistently with all of its moments, not just the first two. When, 
in a simple lifting step, angles drawn from this known distribution are arbitrarily assigned 
to oscillators with different natural frequencies, the fine-scale evolution will move away 
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from the steady state before eventually returning to it. Even though moments are good 
variables for parameterizing the slow manifold, lifting with them is difficult; one needs to 
evolve constrained on the prescribed moments to obtain "mature" distributions, that both 
have these moments and are close to the slow manifold. Observing the evolution of the 
state during the constrained integration sheds light on the nature of difficulty in this lifting 
problem: we observed that during the initial evolution stages correlations between the phase 
angle and the oscillator natural frequency develop. These correlations are destroyed when 
we randomly assign angles from a given distribution to oscillators with different natural 
frequencies. This strongly suggests that different coarse-grained observables that take these 
correlations into account must be sought, so that lifting can be implemented in a simpler 
and less computationally complicated/expensive fashion. We propose that this can be done 
through the framework of the Wiener's polynomial chaos, where both the phase angles and 
the natural frequencies are treated as random variables, and the former is expanded in 
terms of Hermite polynomials of the latter; the expansion coefficients are chosen to be the 
coarse variables [Moon et al, 2005]. 

In our current work we clearly saw that the selection of good coarse-grained observables 
that parametrize the slow manifold, and the construction of a good lifting operator are 
vital for the success and competitiveness of this approach compared to direct simulation. 
For the examples studied in this paper, only a small number of macroscopic observables 
were enough to answer the coarse-grained questions of interest. Equation-free methods 
like the gap-tooth scheme [Gear & Kevrekidis, 2003c] and patch dynamics [Samaey et al, 
2004] can be useful for problems whose dynamics are described in terms of smooth fields 
of coarse-grained observables. The review papers [Kevrekidis et al, 2003; 2004] contain 
detailed descriptions of equation-free methods, as well as references to an extensive set of 
applications where the algorithms have been applied and experience with their properties 
and limitations that has been gained. 
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Figure 1: Qualitative schematic of the fraction of phase- locked oscillators in the Kuramoto model; there 
exist two critical values in K, the coupling strength. Npl is the number of phase-locked oscillators. 
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Figure 2: Bifurcation diagram of the oscillator that desynchronizes first with decreasing K (this is the 
one with the maximum value of \oJi — (wj) |, the " extreme" oscillator). The transition from the full to the 
partial synchronization regime at K = K tp corresponds to a turning point of this oscillator phase angle at 
steady state (0f xp ), where 9t P is the phase angle at the turning point. For the extreme oscillator, a saddle 
state and a stable node state collide in a saddle-node infinite period (sniper) bifurcation at K — K tp . 
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Figure 3: Trajectories of 300 oscillators (K = 0.7) in the phase space of M.2 and M4, started from four 
different sets of values, but with the same natural frequency distribution. Each dot is separated by 0.1 
in time (arrows indicate the direction of increasing time). Initially, g{uS) and f(0) have no correlations. 
In all cases, both moments decrease first (during t < 1), and then gradually increase toward a steady 
state [{M s 2 s ,Mf) = (2.3131xl0~ 2 , 1.6076xl0~ 3 ) marked by a circle]. The latter slow evolution occurs 
along a "slow manifold", which nearly coincides with the curve for Gaussian distributions (dotted line, 
M 4 -3M 2 2 = 0). 
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Figure 4: (Color) Two different random initializations of oscillators (marked by squares) are directly 
guided (shown by dot-dashed black lines) to the slow manifold by solving Eq. (JIJ along with a constraint 
on A^2 for K — 0.7. Each dot is separated in time by 0.1. Poor lifting results from relaxing the constraint 
before the slow manifold is approached (colored lines, constraint lifted at various times t = 0.5, 1, 1.5, and 
6 for red, blue, green, and black line, respectively). Arrows indicate the direction of increasing time. 
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Figure 5: Transient approach to the steady state (K = 0.7) computed through coarse projective integration 
(shown as dots), taking (a) only M. 2 , and (b) both M.2 and M4 as coarse variable(s). Full, fine scale 
integration is shown for comparison (solid line). In the projective integration computation, Eq. Q is 
integrated only for the time intervals marked by dots, and the system is restricted to M.2 (and M4) at 
uniformly distributed times (five dots in each group) . The last three points are used to estimate the future 
value(s); lifting is then performed again, and the procedure repeated. For each lifting step, the same g(u>) 
was used. The inset in (b) shows the blow up near the steady state, where larger projection time steps are 
used for faster convergence. 
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